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ABSTRACT 

This report contain* tha ratult* of a study of computer 
tomographic technique* for reconstruction of three-dimensional 
aerodynamic density fields from inter ferograms recorded from 
several different viewing directions. Kmphasis is on the case in 
which an opaque object such as a test model in a wind tunnel 
obscures significant regions of the inter ferograms (projection 
data). This study was carried out in the Interferometry 
Laboratory of the Department of Mechanical Engineering and 
Applied Mechanics at the University of Michigan during the period 
June 1, 1981 through August 31, 1982. 

The study was carried out in three phases t 

1) The development and empirical study, using simulated 

data from analytical teat fields, of a new method 
called the Iterative Convolution Method (ICM). 

2) An empirical study, using simulated data from analytical 
test fields, of existing methods in which the field is 
repreaet c«d by a series expansion. 

3) The application of computer codes based on each of the 
above methods to analyis of real experimental data in 
the form of aerodynamic inter ferograms provided by 
Georg* Lee of NASA Am** Research Center. 
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List of Symbols 

x,y Cartesian coordinates in the object domain 

r,$ polar coordinates in the object domain 

Cartesian coordinates in Fourier domain 

p,0 modified polar coordinates in Fourier domain, see 

Eq. 3 

8 angle subtended by a ray of 0 th projection and x axis 

in the object domain (measured counter clockwise from 
x axis, 0<8<ir) 

Y Cartesian coordinate y turned counterclockwise about 

angle 8 in the object domain 

Y, 9 coordinates in the projection domain 

r Q radius of the circular region supporting the object 

domain 

p Q radius of cutoff circle in Fourier domain 

f(x,y) true object function (change of the index of 
refraction) 

f(r,<t>) f(x,y) in Polar coordinates 

f(Y,9) projection function, i.e. the Radon transform of 

f(x,y) (order of the interference fringe multiplied 
by wavelength) 

F(£,n) Fourier transform of f(x,y) 

F ( p , 0 ) F(£*n) in modified Polar coordinates 

N number of the available projections 

number of the sampling points (rays) in each 
projection 


M 
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introduction 

The tomographic reconstruction of a 3-D object field can be 
reduced to the reconstruction of its cross sections in a set of 
parallel planes. Therefore, this discussion is concentrated on 
the reconstruction of one 2-D cross-section of a 3-D refractive 
index distribution. Throughout this report the 2-D refractive 
index distribution in this cross-section is denoted as the object 
field and is described by an object function . The object function 
is considered to be space limited with nonsero field only in the 
interior of a circle of finite radius r Q . 

It is also assumed that refractive bending of light rays 
parsing through the object field is negligible. Rays are then 
straight lines and the projection function obtained from the 
inter ferograms is the 2-D Radon transform of the object function. 
This assumption generally is reasonable in wind tunnel diagnos- 
tics, except in the neighborhood of strong shocks. 

The study reported herein is concentrated on the case of 
projection data which are incomplete because opaque aerodynamic 
test models in the wind tunnel block part of the object wave used 
to produce the inter ferogram. Only a few authors, such as Bates 
et. al. 3 , Zien et. al. 4 , and Vest et. al. 5 , have considered 
reconstruction from incomplete data of this type. 

We let n(x,y) be the refractive index distribution of the 
object field and n Q be its value everywhere outside the circle 
of radius r Q then f(x,y) - n(x,y) - n Q is the object function, 

A 

F(£*n) is its Fourier transform and f(Y,0) is its Radon transform, 
i.e., its projection function. In the object domain we use the 
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polar coordinates r,<j> as well as the Cartesian coordinates x,y, 
and in the Fourier domain we use the modified polar coordinates 
p,0, -®£p<*,O<0<ir as well as the Cartesian coordinates g* In 
the projection domain only mutually perpendicular coordinates Y # 0 
(-<x><Yo») are used. If K(Y, 0) is the order of interference fringes 

A 

and X is the wavelength of light then f (Y,0 ) * K(Y, 0 ) X • 

2. The Choice of the Method . 

Before introducing the new method developed under this 
contract for reconstruction from incomplete projections, let us 
briefly review the basic features of standard algorithms for 
reconstruction from complete projections. 

There are three classes of algorithms — Fourier transform 
methods, explicit methods and implicit (series expansion) 
methods . 

The Fourier transform method is based on the Central Section 
Theorem which is summarized by Eq. (1): 

00 ~ 

F(p,0) * / f(Y,6) exp (-i2irpY)dY (1) 

•00 

which indicates that the 2-D Fourier transform F(£,n) - F(-psin0, 
PcosQ ) of the object function f(x,y) along any diameter Is equal 

A 

to the 1-D Fourier transform of the projection function f(Y,0) 
along a line parallel to this diameter. The angle 0 is measured 
counterclockwise from the x axis to the projection lines (light 
rays). Eq. (1) can be derived from the definition of the Fourier 
transform of f(x,y): 
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F(£,n) - / S f(x,y) exp-i2iT( £x+ny)dxdy 

.oo —no 

by substituting 


K * -psinQ 


n - pcos 6 


where 


— °°< p <°o and 


0 <0 <7T 


and then substituting 


Y m -X sin0 + y cos 0 


X » x cos0 + y sin0 


so that we get 


F( p, 0) « / J f (X, Y) exp ( -i 2 tt pY)dXdY 


The quantities p and 0 can be interpreted as modified polar 
coordinates in the Fourier domain. Cq. (1) follows from Eq. (5) 
given that 


f(Y, 0) - / f(X,Y)dX 


because X,Y simply are Cartesion coordinates located parallel and 
perpendicular to the projection direction, respectively. 

A 

Knowing the projection function f(Y,9) we can compute the 
Fourier transform F( p, 9) of the object function f(x,y) along any 
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diameter by using Eq. (1)* Than using ths intsrpolation we can 
obtain F(£, n) from F( p r e) and finally by using ths 2-D invars# 
Fourier transform wa obtain tha objact function f(x,y). 

Explicit mathods ara basad on the invars# 2-D Radon trans- 
form! 


f (r,<j> ) 


1 Y 2 r 9f(Y,0) dYd6 

1? i/2-- rsin( *- S) - v 


(7) 


where, instead of Cartesian coordinates x,y, tha polar coordinate 
r, <J> have bean used. Radon's formula, Eq. (7), directly gives 
tha object field f(r,$) by operating on the 

A 

function f(Y,0 ). The most frequently used approximation 
for evaluating integral (7) was proposed by Ramachandran and 
Lakshminarayanan 6 and results in what is called the convolution 
method . It has been shown that by using the Fo u rier Transform 
(p, 0) we can write 

tt/2 » 

f(r,<f>) * J J | p | F ( p , 0 ) exp[i2Trp rsin(<J> -6 ) ] dpd0 (8) 

-v/2 

instead of Eq. (7). Now, writing formally 


00 

w( Y) - / | p |exp{i2trpY)dp 

—00 


(9) 


using Eq. (1) and letting 

00 

g (r sinU-0),0> - J f (Y, 0) w(r sin( ) - Y)dY (10) 

— 00 
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we can rewrite Eq. (8) as 


ir/2 

f (r, <(>) ■ / g(r sin($-9),0) d6 (11) 

-tt/2 

This means that if we can compute the weight function w(Y), then 

A 

given the projection function f(Y,0 ) we can evaluate the convolu- 
tion (10) and consequently compute the object function f ( r» 4>) 
with help of the integral (11). However, it is obvious that the 
integral (9) cannot be evaluated because it diverges as p-** . 
Computationally, we replace the limits -» and* by -p Q and p 0 , 
where P Q is some large finite number. Therefore the convolution 
method introduces into the reconstruction an error which is 
dependent on the values of the spectrum F(p,9) outside the 
cutoff circle of the radius p Q . The smaller the harmonics outside 
the cutoff circle are, the less is the error in the reconstructed 
object function. However, F(p,0) can never be exactly zero out- 
side any finite cutoff circle because F(p,0 ) is the Fourier 
transform of a spacelimited object function. 

A 

Ip implicit methods we expand f(x,y) or f (Y,0 ) in a series. 
Since the Radon transform is linear, 

f ( Y,0 ) - l Ea-.b.. (Y,0 ) (12) 

if 

f (x,y) » I I a. , b., (x,y) (13) 

i j 

where b^j(x,y) are generating functions, a^j are their coef- 

A 

ficients and b^j their Radon transforms. There are many ways of 
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choosing the generating functions or their Radon transform (see 
the summary in ref. 1 or 5). Having selected the generating 
functions we obtain the coefficients a^j hy solving Eq. (12) 
for a sufficient number of points (Y,6). Knowing a^j we obtain 
the values of the object function from Eq. (13). 

Unlike Fourier transform and explicit methods, implicit 
methods can be applied directly to incomplete projections by 
selecting the points for computation of the coefficients a^ in 
that part of the projection domain where data are available. 
However, the more restricted the region of available projection 
data is, the more ill-conditioned the system of equations for 
computing a^j is. The application of the series expansion method 
to reconstruction from incomplete projections has been studied by 
Vest et. al. 5 . Their method exploits the analytical properties 
of the Radon transform and probably shows the most effective con- 
struction of the series expansion. However, the best choice for 
polynomials to represent the projection function along the radial 
coordinate Y is not obvious. For the reconstruction from complete 
projections the most suitable ortogonal polynomials would have a 
weight function equal to the envelope function 5 . However, for 
reconstruction from projections which are incomplete because of 
opaque objects in the field, Legendre polynomials with the weight 
function one probably would be the best choice. 

For reconstruction from projections which are limited 
angularly or truncated radially, the best choice is probably 
Chebyshev polynomials of the second kind with the weight function 
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i»c.. 


■ sin 



the polynomial* Q n (Y) * U n+ ^ (Y)/(l-Y 2 )*^ 2 


[(n+1) arccos Y]/( 1-Y 2 ) ^ ^ . 

Implicit series expansion methods are often complex to 
program and have relatively long computation times. 

When one wants to use the Fourier transform or convolution 
method for reconstruction from incomplete projections, he must 
somehow extrapolate the projection function (data) into the 
region blocked by the opaque object. This approach has been used 

o 4 

by Bates et. a 1. and Zien et. al. . Zien extrapolated the data 
across the blocked region by using the fact that the zero-order 
moment of the projection function f(Y, 6) with respect to Y has to 
be a constant for all the projections, i.e., any 0. (See ref. 5.) 
There is no more information in ref. 4 about the creation of the 
missing part of the projection function. However, there exist 
additional mathematical criteria for this extrapolation which are 
not. used in Zien's method. Bates suggested the use of a series 
expansion to replace the missing projection data and set the 
object field equal to zero inside the opaque region in order to 
compute the coefficients of this series. There are some short- 
comings of the Bates method. This method introduces con- 
siderable discontinuity into the object function and 
simultaneously brings the complexity of a series expansion 
method . 

Our objective in this research project was to develop a 
reconstruction technique which works well when the projection 
data are incomplete because an opaque aerodynamic test model 
partially obscures the object wavS used to form the holographic 
inter ferogram. Based on our review of available techniques we 
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concluded that the beet technique probably is either a series 
expansion method or some iterative method. The series expansion 
method reported in ref. 5 probably is the best general approach 
because it "fills in" the region in which data are obscured by 
the opaque object with a function that satisfies all mathematical 
conditions for Radon transforms. In a sense it is a 
generalization of Zien's method. 

Iterative processes appeared to be attractive because of 
their success in the related problem of reconstructing objects 
from their autocorrelations. This is a similar problem in that 
one has Fourier transform data that are incomplete because only 
the real parts are known. A similar iterative Fourier transform 
approach to reconstruction from projection data that are 
incomplete because the range of viewing directions is less than 
180 has been successful^. However, for the reconstruction from 
incomplete projections it would be too cumbersome to use revi- 
sions in the Fourier domain because there is no frequency informa- 
tion. Therefore we decided to base our iterative technique on the 

convolution met)’ ^d. The technique is an extension of an approach 

7 

rticently reported by Braga and Vest . This new iterative convolu- 
tion method (ICM) is discussed in the following section. 

3 . The Iteration Convolution Method (ICM) 

The usual, single-step, convolution method cannot be applied 
when part of the projection data are blocked by an opaque 
obstacle. To obtain an accurate reconstruction the ICM 
incorporates the measured projection data and apriori infor- 
mation about the object function through iterative revisions in 
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both the object and projection spaces. 

3.1. Basic Description 

The ICM can be broken into the following operations: 

(a) Initialisation of the missing projection data. (We must 
complete the projection data in order to use the convolution 
method . ) 

(b) Computation of the object field by using the convo- 
lution method. 

t 

(c) Revision in the object domain. (This revision consists 
of zeroing all the values of the computed object function outside 
the circle of radius r Q which is known to contain the object. It 
may also include zeroing either all negative or all positive 
values of the computed object function inside the circle of the 
radius r_ if it is known that the object function is everywhere 
^oiitiVSJiiiS everywhere negative.) 

(d) Computation of a new project function from the revised 
object function by path length integration. 

(e) Revision in the project domain. (This revision 
replaces the computed projection data by the measured projection 
data everywhere that measured data are available.) 

(f) The iteration process is repeated starting at (b) until 
it is terminated by using some suitable termination criterion. 

To describe the method, we introduce a true object function 
f(x,y) defined over the whole object domain including the region 
occupied by the opaque object. A true object function can be any 
function with discontinuities of the first kind whose Fourier 
transform exists and whose Radon transform exists and equals the 
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measured projection data everywhere these data are available. 

One way of studying the properties of the ICM is to study 
the change of the mean square difference E between the true 
object function and the computed object function from one 
iteration to the next. Let F(5,n) and F(p r 0) be the Fourier 
transforms of the true object function f(x,y) in Cartesian and 
modified polar coordinates, respectively. Let h(x,y) be the 

A 

inverse Radon transform of h(Y, 0), which is the revised 

A 

projection function fQf , 0) in the step ( e) of the ICM, and let 
H(£,n) or H(p , 0 ) be the Fourier transform of h(x,y). Further, 
let us denote the quantities and functions computed in the i-th 
iteration by the subscript i. Then the relation between the n-th 
and (n+l)-st iteration c;.n be described as follows: 


00 00 

E n ■ f J | f (x,y ) - f n (x,y) | 2 dxdy * 

—oo —oo 
00 00 

« J J |FU,n) - F n U,n) l 2 d£dn 

—oo —oo 

* / W / !f{ p, 0) - F n ( p, 0 ) | 2 1 p | dpd0 

0 — oo 

7T «» 

* a R J J |F( P, 9) - F n (P, 9)| 2 dPdQ 

0 - 00 

Tf oo A ^ 

* a n J U f " (Y ' 9 ) * V Y ’ 9 > I 2 dYde 

00 * 

1 a n J J !*<*• 6) -h n (Y,9)j 2 dYd 0 

0 —oo 

(revision in the projection domain) 

■n °° - 

* a n / JI F( P' 9) " H n ( P» 9) I ^P * 30 

0 

- f/lFtp, 9) - H n (c, e) | 2 |p| dpde 

u n -oo —CO 
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^ oo CO 

- c 5 - / J|FU,Ti) - H n (c,n)| 2 dCdn 
-. 00—00 

* 3 J J | f (x,y) - h n (x,y)| 2 dxdy 

— 00—00 

^ 00 00 

> 3 ; /|f(x.y> - f n+1 (x,y)| 2 dxdy 

— 00—00 

(cutoff in the object domain) 
a 

v H 17 

5“ E n+1 


where it has been formally put 

oLJ |p<p,e)-p u (p,e) | 2 1 p | dpde 

*" / ¥ riP(p.0)-P u (p.e) l 2 dpde 

o -® 

7T 00 

b , 0 LJ |P(p,6)-H n (p,9) | 2 |p|dpde 

rriF(P, 9 )-H„(P,e) | 2 d P de 
0 -00 U 

and where the 2-D Parseval's relation between the object and its 
Fourier transform has been used and the 1-D Parseval‘s relation 
has been applied to the 1-D projection function and its Fourier 

A 

transform. f^Y^) is the initialized projection function. 

The ICM approaches the true object function if E n > E n + 1 , 
i.e.# if the mean square difference E decreases. The mean square 
difference E R is greater than E n+1 if a n _> b n , as is clear from 
Eq. (14). From Eqs. (15) and (16) we can presume that this 
inequality will be satisfied if the Fourier transform of f(x,y) 
is closer to the transform of h n (x,y) than to the transform of 
f n (x,y), particularly at the higher frequencies. Unfortunately, 

/v 

because the projection function h (Y,0) has a discontinuity of 
the first kind on the boundary separating the domain of the 


(14) 


(15) 


(16) 


11 


measured and computed projections* in each iteration the ICM will 
approach a less and less suitable true object function as its 
higher harmonics beyond some cutoff circle of the radius P Q play 
a more and more important role. However, in practice we restrict 
considerations to the interior of this circle by evaluating the 
integral (9) with finite limits +p Q . (In fact the value p Q is 
determined by the number of sampling points used to represent the 
projection function? see Sec. 3.4.) Therefore, unlike the itera- 
tive Fourier transform method for reconstruction from angularly- 
limited projections the ICM does not converge when recon- 
structing an object field surrounding an opaque obstacle. There 
is always an iteration beyond which the ICM will diverge. 

The divergence of the ICM will occur even when the object 
function is represented by a finite number of sampling points, 
thereby eliminating a real discontinuity of the revised 
projection function. The finer the sampling chosen, the sharper 
the difference between regions of simulated and measured projec- 
tion data will be, so an earlier and faster divergence of the ICM 
can be expected. On the other hand, if a coarser sampling is 
used a poor approximation to the projection function will result 
and consequently a higher reconstruction error will be obtained. 

In conclusion, successful application of the ICM for 
reconstruction from incomplete projections will have a strong 
dependence on two factors : 

(i) Selection of a suitable number of sampling points, and 
(ii) Determination of the optimum point at which to 
terminate the iteration process . 
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Finally, we note that the situation is different for 
reconstruction from angularly limited projections, for which the 
true object function is given ambiguously so that the ICM does 
not have to select the true object function from an infinite set 
of functions. In this case the ICM usually will converge to an 
object function close to the true object function whicfi can be 
recovered from the frequencies inside the cutoff circle of the 
radius P Q - This has been proved for iteration Fourier transform 
method 2# 8, 9 > 


3 . 2 Criteria for Number of Sampling Points 

Evaluation of inter ferograms yields values of projection at 
discrete locations which generally are nonuniformly spaced. The 
integrals (10) and (li) therefore must be replaced by summations. 
Although the sampling points of the projection function generally 
are nonuniformly distributed we will consider here uniform 
spacing A6 an< j Ay between the individual projections and sampling 
points in each projection (including the blocked region). The 
function values of uniformly spaced points can be obtained by 
suitable interpolation of the nonuniformly distributed points. 
In Sec. 3.4 we will see that the convolution method requires the 
uniform spacing. For M uniformly spaced points over each of N 
uniformly spaced projections we have 


Ay . 



(17) 


A0 » 


£ 

N 


(18) 
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where r Q is the radius of the circular domain of the space limited 
object function. 

Rigorous criteria are available for choosing the number of 
projections N and the number of data points M in each projection 
when an object function with band limit P Q is to be reconstructed. 
Using the Whitaker-Shannon sampling theorem and following Klug 
and Crowther 10 , the intervals AY and A0 must obey the fol- 
lowing inequalities: 


ay 


< l 

■ ^ 


(19) 


A0 < 


2 u P 0 V 3 


( 20 ) 


However, in the current investigation the object is space limited 
rather than band limited. If it has substantial Fourier transform 
components outside a circle of radius P Q , aliasing 11 will occur. 
Substituting (17) into (19) and (18) into (20), the following 
criteria for the number of rays M and number of projections N are 
obtained : 


M - 1 > 4r P_ 
— o o 


(21) 


N + 3 > 2^r Q P 0 . (22) 

We see that there is no upper limit on the numbers M and N. 
However, when the ICM is used for the reconstruction from 
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incomplete projections, then from some value of M, as we dis- 
cussed in Sec. 3.1, the larger M is, the greater the discrepancy 
between the computed and measured projection data. Consequently, 
poor reconstruction can be expected when sufficiently large val- 
ues of M are used. 

Combining Eqs. (21) and (22) gives: 



(23) 


The significance of this relation is discussed in subsequent 
sections . 

3 . 3 Initialization 

To apply the convolution method to the reconstruction of an 
object from incomplete projections, one has to fill in the 
missing part of the projection function. Using polynomials for 
this purpose yields a series expansion method. Instead of this, 
we have used a purely numerical method to fill in the incomplete 
projection data thereby initiatizing the ICM. Three possible 
requirements for selecting initial values to replace the missing 
data are: 

(a) The values must represent a smooth continuation of the 
measured projection data within the limitation of the spacing of 
sampling points. 

(b) The maximum absolute value of this artificial data is 
no larger than is necessary to satisfy requirement (a). 

(c) The zero-order moments of all the projections are 
equal . 
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Three initialization schemes, designated by letters A, B and 
C were studied numerically. The first two requirements, (a) and 
(b), were satisfied in all three versions. Version C did not 
satisfy the third requirement. 

Let f(Y^,6 k ) be the value of the projection function at the 
sampling point ( Y i»®fc)» where l£i^M and l£k<N. In all three 
versions of the initialization the differences A L (0^) and Ap(0^) 
are computed in each projection. For the k-th projection, 

A L<V - *< Y L' e k> ~ *< Y L-l' 0 k> 

( 24 ) 

A F^ 0 k^ * f(Yp, e k ) - ^ Y F+l' 9 k^ 

where i*L and i»F (L < F) denote the measured projection values 
closest to the two edges of the blocked region in the projection 
data. The artificially introduced projection data in the blocked 
region are computed from the recursion formulas: 

■ f<*L + u-i>'V + W/* - fi(Y n -V 

( 25 ) 

f< Y F-l- e k> * + W' 1 + 


where 


A(Y C 


V 


2 ^ Y L + U-1 )' 9 k^ “ f ^ Y F-(£ -1)' e k* + 
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- VV> /*»/(WV - 2(1 - 1,,2 ' (26) 

I t - l,2,3...T[t MAX (9 k )/23, (27) 

WV - p ( 0 k> - t-' 9 k>' < J8 > 


and T[ ] denotes a truncated value. Using these relations, the 
simple intialization of the projection data having the first two 
required properties (a) and (b) is made. This describes 
intitialization scheme C. 

A mathematical requirement for a projection function to be a 
Radon transform is that the n-th moment of f(Y,0) with respect to 
Y have circular harmonics of degree less than or equal to n. in 
particular, if n*0 the moment is constant in all projections and 
equal to the integral of f(x,y) over the entire domain^. In 
version B, the zero-order moment of each projection. 


m{© k ) 


r 


A 


f(Y,e k ) dY, 


(29) 


is computed by using the trapezoidal rule. Further, the maximum 
zero moment m M ^ x is determined. Having found the projection 
which has this maximum moment, the previously-computed values of 
the projection function for each of the other projections are 
revised so that each has the same zero-order moment m(6 k ) equal 
to m MAX . In the k-th projection this is done by applying the 
relations 


“ ^ Y l+£'®k^ + D ( Y J,*®fc) 
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(30) 


*NEW* Y F-Ji' 0 k “ ^ Y F-A , 6 k^ + D * Y £' 0 k J 


where 


D(Y r e k ) 


< n *MAX*" m ( 0 k> 1 * 

1+1+S+2+3+3+. . •+(^ MAX ’-1) ' 


(31) 


l MAX ( 0 k) i« given by (28) and integer l goee again through the 
values given by (27). This primarily changes the computed 
projection function in the center of the blocked region and only 
minimally affects the smooth transition of the initialised 
projection function from the measured to the computed domain. 

The version A uses the same revision of the computed 
projection data; however, the maximum zero moment m^^ is 
selected from the projections having minimum i,e * * rom 

the projections that have the smallest blocked region. 

The probability that the initialised projection function is 
the Radon transform of some true object function is obviously 
higher in the case of the versions A and B thar. in the case of 
the version C. 


3.4. Computation of the Object Function 

Having projection data over the whole projection domain we 
can start to compute the object function. 

First recall that the infinite limits of the integral (9) 
must be replaced by sufficiently large finite values -p 0 and p Q , 
After this substitution, Ramachandran and Lakshminarayanan^ 
showed that a suitable weight function w(Y) for use in Eq. (10) 
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may be exactly represented by discrete aample valuea uniformly 
distributed with apacing AY * l/2p Q . They derived 


w(m(Ay)) » 1/4(AY ) 2 for m 0 

» -l/[-rrm( AY) ] 2 for m “ odd 


(32) 


0 for m 


even 


where m is an integer. Note that the Fourier transform of this 
weight function is a triangular function with period 2p Q and not 
the function | p | . 

Now, aaauming that we alao know the projection function at 
aample points with uniform spacing AY - l/2p Q , we can replace 
the convolution integral ( 10 ) by the summation 

. f(m(AY),0.) , f ((m+p) (AY),6. ) 
g(m(AY) , e^) “ 3 " I — K 3 3 ) 

TT D 

p odd * 

However, to evaluate the integral (11) we need to know the values 
of this convolution function not only at the points m(AY), where 
m is integer, but at the points 


y(AY) * r sin^-Q^) 


(34) 


where, for given r and <p , the number y equals 


y * ( BcosQ^-asinO^) /( AY) 


and 0 , 


k( A0) 


(35) 

(36) 


if k ■ 0 , 1 , 2 ... N - l and B * r sin # 

a - r cos <p 


(37) 
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To find these values of the convolution function instead of 
those given by (33) the following linear interpolation is 
commonly usedt 


g(y(AY),e k ) * (m+l-y)g(m( AY), 6^) ♦ (y-m)g( (m+1 ) ( AY) , 6^) . 

(38) 

Finally, we obtain the object function at any point (r,$ ) by 
replacing the integral (11) by the discrete summation 

f(r,*) - ][ g[ y( r , $) ( AY) , 0v] • ( A0) (39) 

k 

In the ICM we used this computational procedure; however, we 
note that different weight functions also have been applied by 
various authors 12, 13# 14 . All these modifications involve some 
suppression of the higher frequencies of the weight function (32) 
in the interval - p Q to p Q in order to reduce noise due to the 
statistical fluctuation in the data. For our purpose the 
Ramachandran and Lakshminarayanan weight function (32) is 
probably the best choice due to the discontinuities between the 
available and computed projection data created in each iteration 
of the ICM. 


3.5. Computation of the Projection Function 

After computing the object function and effecting its 
revision in the object domain, the ICM goes back to tho projec- 
tion domain to make a revision. The computation of (NxM) points 
of the new projection function is made by evaluating the path- 
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length integral (6) (NxM) times. To evaluate each of these 
integrals, we need L samples of the object function, so we must 
carry out the procedure of paragraph 3. (LxNxM) times to find 
the (LxNxM) values of the object function. 

After considering several methods of numerical integration, 
we selected Gaussian quadrature as the most suitable method for 
our purpose. This method requires a small number of sampling 
points for a given accuracy, and they can be distributed 
nonuniformly as is suitable for this application. We assume the 
field has no discontinuities. 

The most time-consuming operation in the convolution method 
is the computation of the number y given by Eq. (35). For the 
computation of (LxNxM) object points we need to compute Nx(LxNxM) 
values of y because each object point is used in each projection. 
This requires N* computations of the sine and cosine functions. 
We have reduced these computations substantially by exploiting 
the fact that the same distribution of sampled object points only 
turned about k(A0), is needed for computation of the individual 
projections . 

Let any object point which is needed for the computation of 
the projection data in the k-th p-ojection have a value of y 
designated as in the i-th projection, 0 < i, k < N • 1. Then 

y ik - r sin(<J> k -i(A0))/(AY) (40) 

where 

<f>k » 4> 0 + k(A0) (41) 
so that 
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» r ain( <fr Q -(i-k) ( A6) )/( AY) 


(42) 


and by using Eqs. (37) 

u ik * C8coa( (i-k) ( A0) )-asin( (i-k) ( A0) )]/( AY) (43) 

where a, B are Cartesian coordinates of (LxM) sampling object 
points for the computation of the O-th projection. Substituting 
l ■ i-k we get 

U ?.+k k " 8cos(H(A0)) - asin(j l\A0)) (44) 

where l - N < i < N -l. Using arrays of (2N-1) terms for 
cosine, sine and p, we have reduced the need for the computation 
of these functions by a factor of (2N-1)/N 2 . 

Gaussian quadrature has been applied with fifteen sampling 
points (L - 15 ) . 

3 ♦ 6 Termination of the Iteration Process 

In Sec. 3.1 it was noted that we expect the ICM method to 
converge to a reconstruction of maximum accuracy and then diverge 
during subsequent iterations. It therefore is crucial to 
determine a reliable criterion for terminating the iterations. 
Attempts to determine this criterion analytically were 
unsuccessful, so we resorted to an empirical investigation. A 
large number of ICM reconstructions of known object functions 
from incomplete projections were made. In each, the RMS error of 
the computed object functions was calculated on the bounding 
circle of radius r Q at the end of each iteration. Similarly, the 
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RMS error of the computed projection function over the region in 
which measured data are available was calculated. By observing 
the values of these RMS errors as a function of the overall 
accuracy of reconstruction, an empirical criterion for 
terminating the iterations was developed. This is described in 
Sec . 4 . 


4. Numerical Investigation of the Iterative Convolution Method 

An empirical investigation of the ICM was carried out 
in order to study its convergence/divergence properties, to 
assess the effects of the number of projections N and the number 
of data points per projection M on the accuracy of 
reconstruction, and to determine an appropriate criterion for 
termination of the iterative process. In addition, various 
schemes for initializations and for revision of the projections 
and reconstructions during itexation were evaluated. The 
accuracy of this method also was compared to that obtained by a 
series expansion method reported in ref. 5. 

4.1 Procedure and Results 

In order to conduct a numerical study, we had to select both 
test object fields (analogous to the gas density distribution in 
actual application) and the shape and size of embedded opaque 
obstacles (analogous to the cross section of the opaque 
aerodynamic test model in actual application). 

The test object functions f(x,y) were selected to be the 
same as those used in ref. 5 in order to make comparisons of 
accuracy. Each object function is zero outside a circle of 
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radius r Q ■ 1. The analytic expressions of the test functions 
used in the study ares 


f(x,y) » exp{l-3n*-a) 2 +y 2 3)tl-<x 2 +y 2 >: 1/2 (45) 


and 


f(x.y> - } ♦ a.Sexpl^SiOiilVi , 

l-(x 2 +y 2 ) l-(x 2 +y 2 ) 

♦ .xpt^[x 2 ^y-o-6) 2 ] , + g ,5exp( -^ A( r°; 6>2 ) ) 

l-Cx 2 +y 2 ) l-(x 2 +y 2 ) 

(46) 


where a is a parameter of the function (45). We used a ■ 0.2944 

and a * 1.0351. For the former and latter values the test 

function (45) is designated as NO lb and NOld , respectively. The 

test function (46) is designated as NQ2 . 

Function NOlb and NOld are two-dimensional Gaussian 

distributions differing from each other in that NOld has a very 

steep gradient close to one boundary. Function N02 has four 

maxima two of which are twice as high as the others. NOld and 

N02 are more complicated and difficult to reconstruct than 

distributions which would be met in subsonic flows near simple 

test models. However, for the investigation of properties of the 

I CM such test functions are very helpful. 

Each of these test functions has been studied with four 

different embedded opaque objects. Two were circular discs of 

radius r„ s 0.3 and r_ -0.6 both centered at the origin and the 
c c 

other two were triangular objects having the sizes and positions 
shown in Fig. 1 where the sketched circle has radius r c * 0.6. 
The larcer triangular is designated by the letters TB and the 
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smaller one is designated by the letters TS. The smaller 
circular disc is designated as C3 and the bigger one as C6. 

In the n-th iteration, the reconstructed test function 
f n ( x i#Yi) is calculated at uniformly spaced discrete points 
(*i # y^) and is compared with the exact values of the test 
function f(x^,y i ). 30 x 30 sampling points of intersections of 

the square net has been evaluated. 648 of these points are 
inside the circle of radius r Q . As in ref. 5, the average percent 
error in the n-th reconstruction is defined as 


1 

I* J 


I J 


l l 

i«lj»l 


f n (x i' y j ) ~ f (x i' y j } 
max jf (x^,yj) | 


. 100 


(47) 


The maximum value as well as any other value is chosen by ex- 
cluding the region occupied by the opaque object. Both the 
average and maximum errors have been used for comparing the 
results from different iterations and different reconstructions. 


4 . 2 Evaluation of the Numerical Results 

The results of the numerical investigation are displayed in 
a series of plots and tables. Both plots and tables show 
dependence of the average and maximum percent error of the 
reconstructed field as a function of various parameters. In all 
the plots, solid and broken lines are used for the average and 
maximum percent error, respectively. Most of the plots are 
organized as five groups of twelve figures. The figures in each 
group have the same figure number and are distinguished by a 
small letter following the number. 
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The twelve figures 2a to are plots of PE (percent errors) 
versus ITR (number of iterations). Three ways of initializing 
the iteration. A, B and C, which were described in paragraph 3.3, 
are considered. The letter R is used to designate the results at 
the iteration specified by the termination criterion discussed in 
paragraph 4.2.4. Figures 3a to 3£ are also plots of PE versus 
INT? however, the parameter is the number of projections N. The 
number of data points M in each projection is fixed. Figures 4a 
to 4 1 are plots PE versus N and M is again fixed. The errors of 
computed points are always taken from the iteration providing the 
minimum average error. Each computed point is marked by the 
number of this iteration. If this iteration differs from that 
indicated by the termination criterion the result of the latter 
iteration is also shown in the plot. The discrete points of the 
plots for which the computation has been made are indicated by 
small circles or small crosses. The circle is used when the RMS 
object error on the edge circle of the object domain has its 
minimum in the same iteration as the RMS projection error over 
the region of available projection data. Otherwise the small 
cross is used. Figures 5a to 51 are plots PE versus ITR, the 
parameter is M and N is always matched to M, using Eq. (23). The 
iteration selected by the termination criterion is again 
designated by R. Figures 6a to 6o are plots PE versus M for N 
matched according to Eq. (23). These plots are made and desig- 
nated in the same way as those of Fig. 4. 

Each letter connected with the figure number determines one 
test object function and one opaque obstacle in agreement with 
the scheme summarized in Table 1. 
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TABLE 1 



opaque test 

obstacle/ function 

NOlb 

NOld 

NO 2 

C3 

a 

b 

c 

C6 

d 

e 

f 

TS 

g 

h 

i 

TB 

j 

k 

1 

none 

m 

n 

o 


For instance, plots of Fig. 6c are plots PE versus M for the 
reconstruction of the test function N02 with the embedded opaque 
obstacle C3. 

In the next paragraphs the numerical results are discussed 
and evaluated. 
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Figure 2d 
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Figure 2h 
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Figure 4d 


56 











ORIGINAL PAC£ IS 
OF POOR QUALITY 




60 





Figure 4i 
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Figure 5e 
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Figure 6k 
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4.2.1 Different Initialisation* 

The plot* PE versus ITR of Pig. 2* to 22, demonstrate how the 
three initialisation schemes A, B and C (»•• par. 3.3) affact tha 
Itaration procaaa. Figures 2a to 2f don't contain tha raaults 
for initialisation B bacausa for a circular opaqua objact tha 
Initialisation A and B ara idantical. 

First of all wa saa from tha plots that tha ICM is not a 
convergent mathod. Generally, tha procass initially approachaa 
tha trua objact function (tha arror is decreasing)* however, 
aftar soma itaration it divargas. For a fixed initialisation* 
tha position of tha itaration providing tha smallast arror is 
shiftad to tha lowar itaration* tha smallar tha opaqua objact and 
tha M simplar N tha tast objact function. In soma casts tha first 
Itaration may provida tha smallast arror. 

From plots of Figs. 2g to 2t wa easily saa that tha worst 
raaults ara provided by tha initialisation B. Comparing the 
results wa have got by using tha initialisations A and C we find 
out that the initialisation A provides as good or better results 
than initialisation C. This is true for all cases except that of 
Fig. 2d. However, in this case the average errors for both 
versions of initialisation are smaller than It (which is 
excellent accuracy for reconstruction from incomplete projections 
when more than 501 of the projection data is missing). The main 
advantage of initialisation A over initialisation B is its 
substantially higher efficiency. When initialisation A is 
applied, a substantially smaller number of iterations is needed 
to obtain the minimum of the PE versus ITR plot. 



In conclusion, initialization A has been selected as the 
best one. All other computations have been made using it. 

4.2.2 Dependence on the Number of Projections 

The plots of PE versus ITR in Pigs. 3a to 35. and the plots 
of PE versus N in Figs. 4a to 4£ show how the results are 
affected by the number of projections N. Evaluation of these 
plots can be summarized as follows: 

For a fixed number of data M in each projection it is suit- 
able to select N in the interval between N * M and N given by Eq. 
(23). The value of N matched to given M by Eq. (23) is indicated 
on the coordinate axis N by the mark "m" in each Fig. 4. \ 

smaller number of projections substantially decreases the accu- 
racy of the computation and a higher value doesn’t improve the 
results. If only one value of N is to be recommended, it would 
be the value of N given by the condition (23). 

4.2.3 Dependence on the Number of Rays in each Projection 
Having roughly defined the optimum ratio between N and M, we 

studied the dependence of the iterative convolution process on 
the number of data M using the corresponding N specified by Eq. 
(23). This means in fact that we studied how the iteration 
process depends on the quantity of complete input data having 
matched the number of projections N with the number of rays in 
each projection M, rather than how it depends on the number of 
rays M alone. Results of this study are presented in plots of ?E 
versus ITR in Figs. 5a to 5 5 and in plots of PE versus M in Figs. 
6a to 6o. 
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For this series of computations the opaque object C6 was a 
circular object with the radius r c » 0.7 instead of r c ■ 0.6. 
Because these are M uniformly spaced sampling points over the 
whole projection, the radius r c was defined by 


M+1-2A 
r c * M^T~ 


(48) 


where 2A is the number of sampling points in the unblocked part 
of the projection. Recall that the edge points at which the 
projection function drops to zero contribute to M as well as to 
2A. The radius r c 0.7 varied less with M than thf radius r c 
~ 0.6, so r Q - 0.7 was used in this study. Note that for different 
values of M the circular opaque obstacles change size only 
slightly; however, the triangular objects change shape as well as 
size . 

The range of the values of M was selected to be 7 to 33. N 
was computed from Eq. (23). Again it is obvious from the plots 
in Figs. 5 that the ICM is not convergent to the true object 
function. At first it approaches the true object function; 
however, after reaching some iteration, which sometimes may even 
be the first one, the error starts to increase monotonicity. 
With an increasing number of input data, i.e. with increasing M, 
and simultaneously with increasing N, the minimum of the plot of 
PE vs. ITR is shifted towards the lower iterations. The value 
change of this minimum can be better followed in the plots of PE 
vs. M in Figs. 6a to 6&. Such behavior means that there is no 
sense to increase the input data over some particular value. For 
the majority of our test models the best results have been 
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obtained for M approximately equal to 20. Only for the test 
object function N02 with circular obstacles was there non- 
negligible improvement of results for M greater than 20. 

4.2.4. Termination of the Iteration Process 

Since the ICM is not convergent, the determination of a 
suitable criterion for termination of the iteration process is 
very important. As stated in Sec. 3.6, the RMS error of the 
reconstructed object function on the edge circle of radius r Q 
(RMSO) and the RMS error of the computed projection function over 
the region where measured data are available (RMSP) provide us 
with information about the behavior of the iteration process. 
Generally, both these RMS errors, i.e. RMSO and RMSP, have their 
minima in different iterations. This can be seen in the plots in 
Figs. 4a to 4^ and Figs. 6a to 6 i. In these plots small circles 
denote points at which both RMS errors were minimum at the same 
iteration and small crosses denote cases in which they were 
minimum at different iterations. 

We considered several possible criteria for termination of 
the iteration process. Finally we selected a very simple one. 
In order to state this termination criterion let us define the 
following notation. In the j-th iteration the j-th object 
function, RMSOj, the (j+l)-s projection function and RMSPj + ^ are 
computed. RMSP^ = 0 because it is the RMS of error of the 
initialized projection function. Furthermore, let the computed 

projection be that which is best matched to the measured data in 
the J-th iteration, i.e. RMSP J + ^ is minimum. Then for J/l the 
resulting object function is taken from the J-th iteration if 
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RMSPj + 2 >. RMSP j and from (J+l) 8t iteration if RMSPj + 2 < RMSPj. 
If J ■ 1 then the result is taken from the first iteration 
if RMSO 2 > RMSO^ or from the second iteration if RMSO 2 K RMSO^. 

The reliability of this criterion is presented in the plots 
of Figs. 2a to 2 I and 5a to 51 where the iteration selected by 
this criterion is indicated by the letter R. We can see that 
this criterion determines the iteration yielding minimum average 
error with reasonably good reliability. 

4.2.5. Alternative Object and Projection Revisions 

Several types of revisions in both the object and projection 
domains were studied numerically. 

The basic revision (BR) in the object plane was to set the 
object function to zero outside the circle of radius r Q in each 
iteration. If it is known from the physics of the problem that 
the object functions is non-negative, a modified revision (R) is 
used, namely, the object function is set to zero for r > r Q and 
at any points in the region r £ r Q where negative values have 
been computed. Table 2 summarizes reconstructions of all test 
object functions using each revision (R and BR). 

In the Table 2 the upper and lower numbers are the maximum 
and average percent errors in the iteration selected by the 
termination criterion. If there are also values in parentheses 
it means that the termination criterion did not yield the 
iteration with minimum average error. These minimum errors then 
appear in parentheses. The test object function (TOF), opaque 
obstacle (OB) and number of rays M are indicated in the first 
three columns of the Table 2. If it is not designated otherwise, 
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V ( : v. 

the number of projection N is matched to M in accordance with the 
condition (23). 

The information in Table 2 indicates that revision R should 
be used whenever it is applicable. 

Revision R in the object domain also has been combined with 
various revisions in the projection domain during the 
reconstuction process. 

We know that each revision in the projection domain 
introduces a discontinuity of the first kind into the revised 
projection function. Such a discontinuity affects the 
reconstruction process very unfavorably (see the discussion in 
the paragraph 3.1) We unsuccessfully tried to alleviate this 
drawback by using only a fraction of the computed change of the 
projection function between two sequential iterations to define 
the new revised projection function. As we expected, this 
technique needed a substantially higher number of iterations to 
attain the best result; however, it didn't lead to a decrease of 
the maximum and average errors as we had hoped. One investigated 
case is demonstrated in Fig. 7, which shows the plots PE vs. ITR 
for the case when only 6% of the computed change of the 
projection function has been used to define the new projection 
function (Note the plots designated as A6.). For comparison, the 
results obtained without reducing the computed change of the 
projection function are presented in plots designated as A. 

We also tried to reduce the unfavorable effect of the 
discontinuity of the projection function by "shrinking" the 
opaque obstacle. In each projection we added the first and last 
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Table 2 - Percent Errors for Different Revisions in the ICM and 
for the Series Expansion Method of Ref. 5 (UChF) . 


:of 

OB 

M 

BR 


R 


R(N=M) 

Rl 


R2 

VChF (N< 

TO lb 

03 

21 

1.5 


1.5 


1.3 

1.5 


1.5 


9.0 




0. 2 


0.2 


0.3 

0.2 


0.2 


2.0 


C6 

21 

5.7 

(3.4) 

5.7 

(3.3) 

5.9 

5.5 

(3.4) 

3.4 






0.9 

(0.8) 

0.9 

(0.8) 

1.4 

3.9 

(0.8) 

0.8 




m . » 

- 

20 

14.4 

(13.7) 

14.4 

(13.7) 

11.6 

12.9 

(12.7) 

14.5 

(13.7) 





1.6 

(1.6) 

1.5 

(1.5) 

1.9 

1.6 

(1.6) 

1.6 

(1.5) 



TB 

2C 

22.2 

(21.1) 

23.0 

(24.4) 

22.9 

21.5 

(20.7) 

23.3 






2.4 

(2.4) 

2.2 

(2.2) 

3.2 

2.2 

(2.2) 

2.1 



TO Id 

C3 

21 

6.1 

(6.2) 

6.1 

(6.2) 

7.2 

6.0 

(6.0) 

6.2 

(6.0) 

9.1 




0.8 

(0.7) 

0. 5 

(0.5) 

0.8 

0.8 

(0.7) 

1.0 

(0.8) 

1.7 


C6 

21 

11.0 


11.3 


12.4 

10.6 

(9.6) 

7.6 

(11.2) 





1.9 


1.2 


2. 1 

2.0 

(1.5) 

2.0 

(1.9) 



TS 

20 

14.9 

(12.7) 

11.0 

(9.1) 

14.2 

8.5 


11.0 

(9.1) 





1.8 

(1.8) 

1.3 

(1.2) 

1.9 

1.2 


1.3 

(1.2) 



TB 

20 

19.2 

(21.1) 

13.9 


17.0 

12.8 


14.7 






2.4 

(2.3) 

1.5 


2. 5 

1.9 


1.6 



;o2 

C3 

27 

19.5 


19. 5 


20. 7 

19.5 


19.5 


10.5 




1.8 


1.7 


2.0 

1.8 


1.8 


2.5 


C6 

27 

27. 5 


26. 3 


27. 3 

26. 9 


36.5 

(30.5) 

17. 5 




3.7 


3.4 


5.2 

3.9 


6.2 

(4.6) 

4.0 


TS 

25 

23.4 


18.7 

(24.3) 

41.0 

17.9 


18.7 

(24.3) 





2. 8 


2.7 

(2.6) 

3.3 

2.5 


2. " 

(2.6) 



TB 

25 

39.4 

(54.4) 

33.6 

(37.6) 

47. 3 

40.2 

(38.5) 

32.8 

(32.7) 





4.9 

(4.6) 

4.2 

(3.9) 

5.3 

3.9 

(3.7) 

4.1 

(3.9) 



TS 

20 

16.1 

(16.2) 

15.0 

(14.9) 

18.2 

19.4 

(18.0) 

15.9 

(15.2) 





2.9 

(2.8) 

2.5 

(2.5) 

3.5 

2.6 

(2.6) 

2.6 

(2.5) 



TB 

20 

22.6 

(23.3) 

20. 3 

(20.4) 

31. 7 

15.6 

(15.0) 

20.3 

(20.1) 





3.8 

(3.5) 

3.2 

(3.2) 

5.0 

3.3 

(3.2) 

3.3 

(3.2) 
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sampling points from the initialized projection data in the 
blocked region to the set of the sampling points of the regions 
where the measured data were available, and considered the 
function values in all these points to be fixed for the iterative 
revision. In this way we in fact numerically shrank the opaque 
obstacle. The results obtained using this modified revision in 
the projection domain are presented in column R1 of Table 2. The 
modification didn't bring any obvious improvement and therefore 
was rejected. 

We made a final attempt to alleviate the influence of the 
discontinuity of the revised projection function on the 
reconstruction: We transformed the problem of reconstruction 

from incomplete projections (in which the true object function is 
not determined unambiguously) to the case when one projection, 
after initialization of the missing projection data, was 
considered to be a complete projection with fixed data over the 
whole reconstruction process. This approach doesn't admit any 
ambiguity of the true object function. The results obtained with 
help of this method are in the column R2 of the Table 2. 
Comparing the results in this column with those in the column R 
we see that they are more often worse than better, so this 
modification of the revision in the projection domain also was 
eliminated . 

In conclusion, the original revisions described in paragraph 
3.1 were found to be the best choice. However, when revision R 
is applicable, it should be used rather than BR. 
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4.2.6. Summary 

From the empirical study reported in this section it is 
clear that the iterative convolution method (ICM) is not 
convergent when an opaque obstacle is present in the field. In 
general as the iterations proceed the reconstruction initially 
converges toward the true object function, reaches a 
reconstruction with minimum error, and then diverges. For very 
simple object functions the number of iterations to achieve 
minimum error is small. It may even occur at the first 
iteration. As the amount of data, specified by the number of 
projections N and the number of data points per projection M, 
increases the curve of percent error (PE) versus number of 
iterations (ITR) becomes more concave and the point of minimum 
error becomes lower and occurs at a smaller number of iterations. 
However, there is “point of diminishing returns" beyond which 
increasing the amount of data no longer increases accuracy of 
reconstruction . 

Several possible types of iterative revision in both the 
object and projection planes were studied, and the best approach 
was selected. Similarly, a rather reliable criterion for 
terminating the iterative process near the reconstruction with 
minimum error was established. Guidelines for selecting an 
optimum number of projections N for a given number of data points 
per projection M were developed. However, in real situations, 
determination of a suitable number of data points M is more 
complicated . 

In aerodynamic interferograms data are usually recorded at 
the centers of fringes whose spacing is nonuniform and whose 
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number varies from projection to projection. The convolution 
algorithm requires equally-spaced data points in each projection, 
hence the M equally-spaced must be obtained by interpolation of 
the experimental data prior to running the code. In addition, 
the low quality of aerodynamic inter ferograms encountered in 
practice may require extrapolation of projections into the free 
stream. Considerable operator judgment may be required in these 
steps . 

The results of our empirical study verified the obvious fact 
that the number of projections N required for accurate 
reconstruction increases with increasing complexity of the object 
field. The optimum numbers of projections in this study were in 
the range 10 to 30, i.e., views would be taken in intervals of 
6° to 18°. 

Some understanding of the accuracy of the ICM can be obtain- 
ed by studying columns BR and R in Table 2. Clearly it depends 
on the complexity of the object function. Average (maximum) per- 
cent error varied from 0.2 (1.5) for test object field 
N01b + C3 to 4.9 (39.4) for test object N02 + TB. The latter is 
•a worst case and could be improved to 3.2 (20.3) by using 
M * 20, rather 25, and using revision R, rather than BR. 

Column VChF of Table 2 contains results of reconstruction 
using a series expansion method presented in ref. 5. Comparing 
these with the ICM results in Column R, it is seen that the ICM 
generally gives smaller average errors. Lower maximum errors were 
produced by the series-expansion method for object function N<32, 
but the ICM produced lower maximum errors for functions N01b and 


N01d. 


Finally, Fig. 8 indicate* the dependence of average and 
maximum reconstruction errors, attained with the ICM, on the size 
of the opaque object specified by the ratio r c /r Q . 

4.3. Increasing the Number of Projections by Interpolation 

In paragraph 4.2.2 we have recommended how to select the 
number of projections N for a given M. N should be roughly given 
by condition (23) even though for some object functions it is too 
stringent a requirement. However, in no case should the number of 
projection N be smaller than M. 

In some practical applications the number of available 
equally spaced projections may not be high enough. To overcome 
this trouble we tried to increase the number of available 
projections by interpolation. Since the projection function 
f ( Y, 6) is periodic in 6 with period 2 tt (considering Y>0 h«re) we 
assumed that triangular interpolation using harmonic analysis can 
be used. Of course, the periodicity of the projection function 
f(Y ,0) is helpful for interpolation in the 0 direction only for 
the values Y which are outside the radius of the region of the 
blocked projection data. If a line Y <=* const, is broken by the 
blocked region, trigonometric interpolation implies that the 

A 

function f(Y » const, 0 ) is periodic with the period defined by 
each individual interval of available data. Obviously, this im- 
pairs the accuracy of the interpolation. Only a circular opaque 
obstacle centered at the origin doesn't break any line Y * const. 

The usefulness of azimuthal interpolation is indicated by 
Table 3. The first of three columns of this table indicates the 
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the same notation as Table 2. The columns designated by N show 
the number of available "measured" projections. The next column 
presents maximum and average percent errors for reconstruction 
using these "measured" projections. The change in errors obtained 
by increasing the number of projections by a factor of 2 or 3 
using triangular interpolation can be seen by noting the errors 
in the columns designated by 2 or 3. In all cases the ICM was 
used with revision R in the object domain. 

From the results in the Table 3 it is obvious that the 
increase of an unsatisfactory number of projections by 
interpolation can be very useful, particularly when a centered 
circular obstacle blocks a part of the projection data. 

5 . Program ICOM; The Iterative Convolution Code 

ICOM is the final computer code for the iterative 
convolution method described in Sec. 3. The source program was 
written in two versions of FORTRAN IV. One version is for 
execution on a Digital Equipment Corporation LSI 11/23 
microcomputer. Execution on this computer required from several 
minutes to one or two hours, depending on the amount of data and 
complexity of the field. The second IBM version was in stan- 
dard FORTRAN IV and was executed under MTS (Michigan Terminal 
System) on an AMDAHL 470V/7 in less than 16 seconds of cpu time 
in a 1 1 instances. Both versions of the code ICOM are included in 
the appendix of this report. 

The two-dimensional object function f(x,y) is always 
considered to be zero everywhere outside the circle of radius 
r Q . This function .is defined in a two-dimensiona 1 right-handed 
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Cartesian coordinate system x,y with origin at the center or the 
circle. Modified polar coordinates R,$ also are used in the 
object domain. These coordinates are related to the Cartesian 
coordinate by the relations 

x = R cos <J) y = R sin <{> (49) 

and are defined in the interval - 00 <^R <_ a0 and 0 4> <_ tt . The 

/v 

individual projections f(Y,8 =constant) are described using a 
coordinate system Y which is perpendicular to the optical rays 
used to form it. The first projection (6 = 3) is formed by 
optical rays parallel to the x axis; the coordinate Y is then 
parallel to the y axis. Every other projection is described by 
rotating the axis Y counter-clockwise with respect to the xy 
coordinate system. 

The detail ed input directions and the output description 
of the program ICOM follow. 

5.1. Input 

The device name for the input is the number 5. The whole 
input consists of 3 + 2N records. The first record is an 
arbitrary text of up to 120 characters. If the IBM version is 
used the first character must be a blank. 

The second record comprises the ten following numbers: 
M PLOT, MSGNV, IC, N, M, NET, NET2, D., WL., GDC. where integers 
and real numbers are without and with the periods, respectively. 

The integer MPLOT determines the optional output. 

MPLOT >0 In this case the computed values of the object 
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function are arranged into a plot of the investi- 
gated cross-section of the object field. The 
values of the object function are computed and 
printed at the intersections of a square net. 

The vertical and horizontal lines of this net 
are parallel to the y axis and x axis, respec- 
tively. 

MPLOT ■ 0 In this case the values of the object function 
are computed at the same points as in the pre- 
vious case; however, the computed values are 
arranged into a table. 

MPLOT < 0 The values of the object function are computed 

at the intersections of equally angularly spaced 
diameters of the circular object field with 
equally spaced concentric circles. The computed 
values are arranged into a table. 

For more detailed information about the output arrangement 
see Sec. 5.2 

The value of the interger MSGNV determines which of two 
possible revisions in the object domain is used. 

MSGNV » 0 This alternative is prefered when it is known 
that the object function does not change sign 
anywhere in the object domain. The ICOM uses 
revision R (See Sec. 4.2.5.). 

Otherwise the object function can change its sign in the object 
domain. Then ICOM uses the revision BR (See Sec. 4.2.5.). 


to increase the number of projections by the trigonometric 
interpolation discussed in the paragraph 4.3. The trigonometric 
interpolation is not a part of the program ICOM. It has been 
coded separately in the program INTEP. The integer IC is selected 
differently to be nonzero only if we use both the programs ICOM 
and INTEP. The integer IC is set equal to a positive multiplica- 
tive factor by which the number of projections is to be increased. 

The integer N equals the number of available projections. 
They must be distributed with equal angular intervals. The 
maximum allowed value of this integer is 43. 

The integer M equals the number of equally spaced sampling 
points over each entire projection including the edge points on 
the circle of radius r Q and the sampling points inside the 
blocked region. The maximum allowed M is 30. 

The integer NET is connected with the output. It specifies 
the density of the sampling points at which the resulting object 
function is computed. If MPLOT 0 and the reconstructed object 
function is computed at the intersections of the square net 
overlapping the investigated cross-section of the object field 
then NET x NET is the number of the sampling points, i.e., inter- 
sections of the square net, at which the object function is 
computed. The side of the square net equals the diameter 2r Q of 
the circular region of nonzero object field. If MPLOT < 0 and the 
reconstructed object function is computed at the intersections of 
the "polar net", then NET defines the number of the sampling 
points, i.e., intersections, on each diameter. The maximum al- 
lowed NET is 60. 

The integer NET2 completes the information about the 


distribution of the output sampling points when MPLOT < 0. In 

this case NET2 defines the number of the equally angularly 
spaced diameters of the output "polar net". The maximum 
permitted number of these diameters is 60. In connection with 
MPLOT > 0 the value of the integer NET 2 is arbitrary and has no 
influence on the computation. However, some value must be 
assigned to NET2 in the input record. 

The real number D. is equal to the diameter 2r Q of the 
circular region containing the nonzero object function and may be 
expressed in arbitrary length units. 

The real number WL- is the wavelength of light in the same 
units as the diameter D if the goal of the computation is the 
reconstruction of an object function from measured projections. 
However, for testing purposes we can set WL. * 0. to indicate 
that we wish to test the accuracy of reconstruction of a test 
object function TF specified as a function subprogram of the main 
program ICOM. In this option the percent errors of the 
reconstructed object function, rather than the function itself 
are computed . 

The final real number GDC. of the first record enables one 
to compute the object function f(x,y) as the relative variation 
of the index of refraction n(x,y) 

f ( x, y) * n(x,y) - n Q (50) 

or as the relative variation of the density p(x,y) 


1 
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f(x,y) 


(51) 


p(x,y) - p Q 

The former representation is obtained by setting GDC. equal to 
zero and the latter is obtained by setting GDC. equal to the Glad- 
stone-Dale constant. Then the computed density is in units of 
the Gladstone-Da le constant to the minus one. Note, that the 
selection of the constant WL. overrides the selection of the 
constant GDC. 

The third record consists of 2N integers A^, A 2 , •••^ 2 ^ which 
define the regions of the individual projections that are blocked 
by an opaque obstacle. Generally, each projection consists of 
three regions - a central region of blocked data and two regions 
of available data, one to either side of the opaque object. The 
blocked region in the k-th projection is defined by the pair of 
integer variables A 2 ^_^» A 2 ^ where k * 1, 2,...N. For any pro- 
jection the odd variable, &2X-1' specifies the number of data 
points lying to the left of the opaque object (i.e. smaller Y) 
and the even variable, specifies the number of data points 

lying to the right of the opaque object (i.e. larger Y). In the 
first projection Y * y and each next projection is rotated 
counterclockwise in equal angular intervals through 180®. 

The boundary of the blocked region generally lies between 
two sampling points. We have found that the best accuracy is 
obtained if the data are smoothly extrapolated to the first 
sampling point inside the blocked region. This point then is 
included with the measured data. 

Because the samp! ing points at the edge of the field (r »r ) 


111 


must be included in the data set no integer A n can equal zero. 
An error message will result if A n ■ 0. Another restriction is 
the A2K-1 + A 2k - M+l. This sum equals M+l if no opaque 
object is present. In this latter case one central sampling 
point should be included in both the "left-" and "right-hand" 
regions of data. 

2N 

The remaining 2N records are the £ A n real values of the 
measured projection function expressed as interference fringe 
orders. It sould be reemphasized that the number of values of 
the projection function must equal the number of sampling points 
in each projection specified by the pair of variables A n . 

Each of next 2N records is related to one value A n of the 
third record. The fourth record consists of Aj^ function values 
of the "left hand" region of the first projection, the fifth 
record consists of A 2 function values of the "right hand" region 
of the first projection and, generally, the (2k + 2) or the (2k + 
3)-rd record consists of or A 2 ^ function values of the 

"left hand" and "right hand" regions of the k-th projection, 
respectively. The individual function values of each record are 
always ordered with increasing values of the coordinate Y, i.e. 
from the edge value of the "left hand" region towards the edge 
value of the "right hand" region. Therefore the first value of 
each record for which A n has an odd subscript and the last value 
for which A n has an even subscript must be zero. (The function 
values at the edge points are zero.) 

Obviously, these 2N records are not required when in the 
second record the option V7L. * 0. is chosen. 
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5.2 Output 

The device name for the output is the number 6. The 
arrangement of the output supposes that a printer with at least 
128 printable columns is available. Before printing the 
reconstructed object function at selected sampling points, the 
following data are printed: The information in parentheses is 

printed optionally and the printed numbers are replaced by 
symbols if they have been introduced or simply by X in the output 
format : 


RECONSTRUCTION FROM INCOMPLETE PROJECTIONS 
(BY THE ITERATION CONVOLUTION METHOD) 

An arbitrary text up to 120 characters 


ZERO OBJECT FIELD OUTSIDE THE RADIUS: 

r o 




INTERPOLATION COEFFICIENT: 

IC 




(WAVELENGTH : 

WL) 

only 

for 

WL f 0 

(GLADSTONE -DALE CONSTANT: 

GDC) 

only 

for 

GDC f 0 


(NO SIGN CHANGES OF THE OBJECT FUNCTION ARE PERMITTED) only for 

MSGNV - 0. 

NUMBER OF PROJECTIONS : N NUMBER OF RAYS IN EACH PROJECTION: M 

(PROJECTION, PROJECTION FUNCTION IN THE ORDERS OF INTERFERENCE 
FRINGES: 

1 x«#x«... 

2 

• • 

• • 

. . ) only for WL ^ 0 

ITERATION, RMS OF THE BOUNDARY OBJECT ERROR AND THE AVAILABLE 
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PROJECTION ERROR: 


x 

X 


X 

X 


1 
2 
• 

• 

RESULT IS FROM ITERATION: X 

MISSING PART OF THE PROJECTION DATA IN PERCENTS: X 

(NUMBER OF OBJECT POINTS, MAXIMUM AND AVERAGE ABSOLUTE ERRORS IN 
PERCENTS : 

x x x ) only for WL * 0 

MAXIMUM ABSOLUTE VALUE OF THE OBJECT FUNCTION: X 

(X and Y COORDINATE MULTIPLIED BY : X) only for MPLOT 1 0 
(THE RADIUM R IS MULTIPLIED BY: X) only for MPLOT < 0 

(THE OBJECT FUNCTION MULTIPLIED BY: X) only for WL 0 

After these data the print of the values of the 
reconstructed object function in the refractive index 
representation (WL j 1 0, GDC * 0) or in the density representation 
(WL f 3, GDC * 0) or the print of the percent error of the 
reconstructed object function (WL * 0) follows. The arrangement 
of this part of the output data depends on the value of MPLOT. 

Let us first suppose MPLOT > 0. Then the values of the 
reconstructed object function or its percent errors are computed 
and printed at the NET x NET intersections of a square net whose 
vertical and horizontal lines are parallel to the y- and x-axis, 
respectively. The width of this net equals the diameter 2r Q . 
The values of the x- and y-coordinates are printed, with an 
accuracy of 1% of the radius r Q , on the upper horizontal and left 
vertical margin of the plot, respectively. The values of the 
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reconstructed object function or its percent errors are also 
printed, rounded to 1% accuracy of the maximum value of the 
reconstructed object function. Four asterisks are printed at the 
sampling "points" in the region occupied by an opaque obstacle 
and also outside the circular region of radius r Q (including the 
edge points). We know that the object function is zero on this 
circle and everywhere outside it. Using this printing technique 
we obtain not only the values of the reconstructed object 
function at the sampling points but also the real geometrical 
image of the circular croas-section of nonzero object field with 
an embedded opaque obstacle. The geometrical image of the real 
cross-section can be formed only up to NET ■ 30. For NET > 30 
the width of the printer is not sufficient for making this 
geometrical image. In this case the output is arranged in such a 
way that the geometrical image or the investigated cross-section 
would be obtained if each even line of the plot were translated 
to the right and attached to the end of the preceding odd line. 

The function values at the sampling points of the square net 
also can be arranged into a table. The x and y coordinates are 
indicated at the top of this table and the values of the 
reconstructed object function, or its percent errors, are without 
indication. This kind of the output is obtained by setting MPLOT 
- 0 . 

By selecting MPLOT < 0 a table similar to that described 
above is printed, but the values of the reconstructed object 
function or its percent errors are computed at the sample points 
of the "polar net." In the top of the table the designation DEG 
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and R are instead of X and Y. The NET printed values for DEG ■ 0 
( $ ■ 0) are those along the x axis. The polar axis is then 
rotated counterclockwise in increments of 180/NET2 degrees. 

Different examples of the output are in the appendix to this 
report. 

6. Program INTEP 

In Sec. 4.3. it was shown that it can be useful to increase 
number of available projections by the trigonometric interpola- 
tion in the azimuthal direction. The program INTEP is used for 
this purpose. 

The source program has been written in FORTRAN IV in both 
modified version as the program ICOM. Copies are enclosed with 
this report. 

The input device is designated by the number 4 and the ouput 
device by the number 5 (which is the input device name of program 
ICOM). The input of the program INTEP is exactly the same as the 
input of the program ICOM. Only the quantity IC of the first 
record is set equal to 1, 2, 3... instead of 0 if we want to 

increase the number of the available projections 1, 2, 3. ..times, 

respectively. The output of the program INTEP is arranged in such 
a way that it can directly serve as input to the program ICOM. A 
file containing the output of INTEP can be used directly as an 
input file for ICOM. The number of projections printed in the 
output of the program ICOM is the number of projections after 
interpolation. The number of projections before interpolation is 
obtained by dividing this number by the printed value of IC 4 0. 
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Note that the program INTEP also contains the function 
subprogram TF. For the choice WL » 0 this subprogram must be the 
same for both the main programs ICOM and INTEP. 

7 . ICOM Reconstruction of Data from Aerodynamic Experiments 

The codes developed during this investigation were used to 
reconstruct aerodynamic density fields using data supplied by 
George Lee of NASA Ames Research Center. We were supplied with 
photographs of dual-plate holographic interferograms of subsonic 
windtunnel flows past a circular cylinder of radius r fl » 25.4 mm 
with a hemispherical nose. 

One data set was for the axisymmetric flow when the cylinder 
was set at zero angle of attack. The second set was for the 
asymmetric flow when the cylinder was set at an angle of attack 
a - 5.5° . 

The density distribution in the flow fields were described 
in our analysis using a Cartesian coordinate system with the z 
axis parallel to the axis of the wind tunnel and pointing in M»e 
direction of flow. The origin was placed at the stagnation point 
of the cylinder at zero angle of attack and the density distribu- 
tion was reconstructed in several planes z ■ constant. Similarly, 

for the asymmetric flow about the cylinder at angle of attack a ■ 

0 

5.5 , the origin of the coordinate system was located where the 
hemispherical nose had a single point of tar.gency to the xy plane 
normal to the flow direction. Thus for z > o the test model 
blocked part of each projection. 
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7 . 1 Axisymmetric Flow (a ■ 0) 

The iterative convolution method (ICM) was developed for 
general computer tomographic reconetruct ion of asymmetric flow 
fields from sets of interferograms when an opaque object blocks 
part of the field. However, it can also be used for the 
reconstruction of axisymmetric flows from a single inter ferogram. 
We first tested the program ICOM by using it to reconstruct such 
an axisymmetric flow field and comparing the results with a 
reconstruction from the same data computed with a code specifical- 
ly based on the Abel transform for reconstruction of axisymmet- 
ric fields 

The interferogram, recorded in an air flow at a Mach number 
M - 0.6, is shown in the appendix. In theory, the interferogram 
should be perfectly symmetric about the z axis and the free 
stream should contain no fringes. However, the actual 
interferogram is not exactly symmetric, and the free stream 
region contains many fringes, which can be considered to be noise 
in the data. The extraneous fringes may be caused by turbulence, 
deformation of the windows or tunnel structure, and by emulsion 
shrinkage or swelling in the two holograms used to form the 
interferogram. It therefore was necessary to use judgement and 
information from a numerical prediction of the flow to: (i) 

determine the zero-order fringe, (ii) define the radius r Q beyond 
which the density change can be assumed to be zero, and (iii) 
decide which half of the interferogram is more accurate. 

Fig. 9 shows the input data which were used in the form of a 
plot of fringe order versus radial coordinate Y. The solid 
curves are interpolations among fringe orders and the dashed 
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curves are smooth extrapolations into the free stream. Each 
curve is labled with the z location of plane of reconstruction in 
millimeters. Short vertical lines denote the edge of the opaque 
test model when it blocks part of the field of view. 

Figs. 10 a, b and c show the reconstructed density distribu- 
tion (FUNCTION) as a function of radius (POSITION) in centimeters 
in the planes z * 0, 12.7 and 25.4 mm respectively. The dark 
curves were produced by the code ICOM and the light curves were 
produced by a special code for reconstruction of axisymmetric 
fields. The agreement between the two is excellent. 

O 

7 . 2 Asymmetric Flow Z. - 1 1 i. 

We were supplied with inter ferograms of this flow for M = 
0.6 and M * 0.8 recorded with viewing directions every 10° over 
the full 180° range. Because this flow is symmetric about the 
plane y * 0, this set of data should be redundant, i.e. the 
inter ferogram for 6, 0^9^90°, should be identical to, or a 

o 

mirror image of* that for (183 - 9). 

These sets of data contained considerable noise due to a 
variety of experimental problems. Density field reconstructions 
could be made only after a number of ad hoc assumptions were made 
to interpret the inter ferograms : 

(i) Because the zero order fringe was not identifiable, we 
postulated that the fringe order at the stagnation point x = y = 
z ~ 0 was the same in each interferogram at a given fiach number. 

(ii) We further postulated that the order of the fringe at 
x = y = z = 0 was 3 for M = 3.6 as i.n the case of * 3, M = 0.6. 
The corresponding fringe order was taken to be 5 for M = 3.8. 
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This value was arrived at by estimating the zero order fringe in 
each inter ferogram to be that which is the first to greatly 
broaden and diffuse into the free stream. The orders estimated 
in this way were averaged to arrive at the value 5. 

(iii) The order of interference fringes was assumed to 
decrease monotonically from that atx=y=z=0. 

Because the ICM performed well in the empirical studies of 
noise-free data described earlier in this report, we believe that 
inaccuracies and artificial fluctuations in the reconstructions 
of the experimental data are due to experimental noise and our 
consequent inability to evaluate the inter ferograms with good 
precision. In order to further evaluate the ICM we reconstructed 
the same data sets using a series expansion algorithm. These 
results are described in Sec. 10. 

Having calibrated all the inter ferograms we had to find the 
projection of the axis z into each inter ferogram. In each 
inter ferogram the projection of the mechanical axis of the 
cylinder could be determined with the help of the shadow image of 
the cylinder. Knowing this projection and using the simple 
relation 


tan a = tan a sin 0 , (52) 

O 

where a is the angle a = 5.5 projected into the 0-projection, we 
have constructed the projection of the z axis in each interfero- 
gram. The corresponding line had tc go through the tangential 
point of the plane z = 0 and the shadow image of the cylinder, 
i.e. point x = y = z = 0. 


124 


We reconstructed the density field in the cross-sections z * 
0, 1.27 and 2.54 cm. The dependence of the interference order 10 
of the projection function on the position Y is shown in the sets 
of plots. For the flow at M = 0.6 these plots are in Fig. 11a to 
1 1 j and for the flow with M = 0.8 they are in Fig. 12a to 1 2 j . 

o o o 

The designations a, b, c, d....j correspond to 9 = 0 , 10 , 20 , 
30°(-30°)...90° (-90° ), respectively. In these plots the circles 
and crosses represent the measured values and the curves show our 
visual interpolation and extrapolation. For z = 0 there are only 
positive values of interference order and projection data are 
available over the whole projection domain. For z * 1.27 and 
2.54 cm the orders of the interference fringes are only negative 
and the central part of the projection function is blocked by the 
test model. The projection functions for z = 2.54 cm always lie 
below those for z * 1.27 cm. For interpolation we tried to 
construct symmetric curves for 9 * 0° • For 2 = 0, following 
the physics and mathematics of the problem we tried to keep the 
area under the curve 10 vs. Y roughly constant for all 
projections. In each cross sectional plane, z = 0, 1.27'and 2.54 
cm, the radius beyond which the density difference is essentially 
zero was estimated by visual inspection of the inter ferograms and 
extrapolated data. 

Equally distributed values of the projection function taken 
from the interpolation/exterpolation curves of Figs. 11 and 12 
were used as input data for the programs INTEP and ICOM. We used 
the program INTEP to double the number of available projections. 
(The program INTEP was not used for the reconstruction of the 
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object field in the cross-section z ■ 0 and M ■ 0.8.) The 

printed input data and the printed results are in the Appendix 

YT " • c • ' ' 

pages A39 to A 106. Each r -construction is presented two times. 
One presentation is with the option MPLOT > 0 and the other for 
the option MPLOT < 0. (See Sec. 5.2.) The former alternative 
gives a map of the reconstructed density change Ap(x,y) ■ p(x,y) 
- P Q in the investigated cross-section. The latter option tabu- 
lates the reconstructed density difference Ap (R, 0 ) = p(R, 0) - p 

o 

in polar coordinates with angular step A8 * 5°. In both cases 
the constant P Q equals the density in the free stream. For 
better illustration the sets of plots in two series of Figs. 13a 
to 13g and 14a to 14g show the graphical presentation of these 
results for the Mach numbers M * 0.6 and M * 0.8, respectively. 
Figures designated by letters a, b, c, d, e, f and g are plots 

o o o o o o o 

Ap vs.R for <t> = 0 ,15 , 30 , 45 , 60 , 75 and 90 , respectively. 
Only the solid lines have been reconstructed by the ICM. The 
broken lines are reconstructions from the same projection data by 
a series expansion method which will be discussed later in Sec. 
10. The three solid curves correspond the planes z = 0, 1.27 and 

2. 54cm. 

8 . Program SEEX 

In table 2 of Sec. 4 reconstructions obtained by the 
iterative convolution method (ICM) were compared with 
reconstructions, based on the same data, obtained by application 
of a series expansion method described in ref. 5. In this method 
the unknown object field is represented by a series of ' 'gendre 

polynomials in the radial direction and by trigo; ometric 

PRECEDING PAGE BLANK NOT FILMED 
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function# in the azimuthal direction. A new computer code 
suitable for analysis of the type of incomplete data considered 
in this report was written. 

The source code for SEEX is included in the appendix of this 
report. It is written in FORTRAN IV. A version for the LSI 
11/23 microcomputer was not prepared because storage requirements 
are too great. The coordinate system, object function and 
projection functions are defined the same way for SEEX as for 
ICOM. (See Sec. 5 and Eqs. (50) and (51).) 

8.1 Input 

The designation of the input device is 5. The input 
consists of seven records. 

The first record is an arbitrary text of up to 120 
characters; the first character must be a blank. 

The second record is read using a NAMELIST specification 
statement called NAM and contains the values of the following 
eleven variables: D, WL, GDC, MSCNV, EPSL, NN, N 1 , M*, MPLOT, 

NET and NET2. Except for EPSL, NN, N 1 and m' all these variables 
have the same meaning as in the program ICOM which is described 
in Sec. 5.1. 

The real value EPSL specifies a relative tolerance for 
determination of the rank of the matrix of the subroutine DLLSQ. 
In the main program this subroutine solves the system of linear 
equation of a least squares method 16 * 

The integer NN equals the number of available projections. 
Generally, these projections are arbitrarily distributed. The 


maximum permissible value of NN i* 36. 

Tho integers N and M respectively define the numbers of 
the Legendre polynomials and the number of harmonics used to 
approximate tho projection function. The maximum permitted value 
of N or M is 12. 

The third record consists of NN angular values of the 
individual projections in degrees. Note that the angle is 
subtended by a ray of a projection and the x axis in the object 
domain, i.e. by the projection plane and y axis. This angle is 
always measured counterclockwise from the x or y axis to the ray 
or the projection plane, respectively. The values of are real 
numbers which can be arranged arbitrarily. However, their order 
defines the order of the individual projections. All subsequent 
records must be consistent with this order. 

The fourth record consists of NN integers defining the total 
number of the sampling points in each individual projection. 
They have to be ordered in agreement with the order of the 
individual projections defined in the third record. The maximum 
number of the samplinq points In one projection is 25. The 0 
points on the circle of radius r Q which defines the object domain 
may have to be among the sampling points. 

Tho fifth record consists of 2(NN) real numbers equal to the 
linear coordinates Y of the boundary of the Mocked region in the 
individual projections. Two of these values Y * Y B(2k-l) 
Y n ( 2k ) ' k * 1, 2...NN belong to each projection. The 

individual couples of these values have to be ordered in 
agreement with the order of the individual projections in the 
third and fourth records. For each projection the smaller value 
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i. 


Y q has to be put first. Therefore, if the opaque object overlaps 
the origin of the coordinate system x,y the first boundary value 
Y q is always negative and the second one is positive. If no 
opaque object is embedded then the fifth record consists of 2(NN) 
zero values. 

The sixth record contains the coordinates Y (real numbers) 
of the sampling points in the individual projections. The order 
of these coordinates in the individual projections can be 
arbitrary but the order of the groups of these coordinates 
belonging to the individual projections must be the same as that 
of the individual projections in the previous records. The 
number of the coordinates in each group is defined and must be 
equal to the corresponding integer of the fourth record. 

The seventh record has the same number of real values as the 
sixth record. These are the values of the projection function 
expressed as orders of interference fringes. They must be 
arranged in exactly the same way as the corresponding coordinates 
of the sixth record. The seventh record is optional. When WL of 
the second record equals zero there is no need for the seventh 
record . 


8 . 2 Output 

The designation of the output device is 6. The output of 
program SEEX is organized like the output of program ICOM. There 
are only the following changes between these two outputs: 

1) The printed title is now: 

RECONSTRUCTION FROM INCOMPLETE PROJECTIONS 
(BY THE SERIES EXPANSION METHOD) 
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2) There is no interpolation coefficient IC in the program 


SEEX; therefore no information about it is printed. 


3) After the optional printing of the G1 adstone-Da le 


constant the following is printed: 


TOLERANCE PARAMETER: 

NUMBER OF PROJECTIONS: 

TOTAL NUMBER OF SAMPLING POINTS: 
of the third input record) 


EPSL 


NN 



(where I k are integers 


ORDERS OF LEGENDRE POLYNOMIALS AND TRIGONOMETRIC FUNCTIONS IN 
SERIES: N' M' 

PROJECTION, ANGLE, NUMBER OF SAMPLING POINTS, BOUNDARIES OF 
OPAQUE OBSTACLE, RELATIVE POSITIONS OF INDIVIDUAL SAMPLING POINTS 
(PROJECTION FUNCTION AT THESE POINTS IN ORDERS OF INTERFERENCE 
FRINGES) only for WL ^ 0 


k 6 k X k Y B ( 2k- 1) Y B(2k) 

(2Y lk /D) (2Y 2k /D). . . (2Yj k /D) (normalized 

k 

coordinates of the sampling points of k-th projection) 

A A /V 

( f (Y^# 6 k > f(Y 2k' 6 k ) - * ,f{Y I k' 9 k }) onlY for 

k 

WL / 0 


this is repeated for all k * 1,2, 3... NN 

4) Then the same print as in the program ICOM starting at 
the optional row 

NUMBER OF OBJECT POINTS, MAXIMUM AND AVERAGE ABSOLUTE ERRORS IN 
PERCENTS: is used. (See Sec. 5.2.) 
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The examples of different outputs as well as inputs can be 


found in the Appendix pp. A35 to A38 and A107 to A154. 


9 . The Reconstruction of Test Fields by SEEX 

The program SEEX was used to reconstruct test object fields 
NOlb, NOld and N02 in combinations with opaque obstacles C3, C6, 
TS and TB. In order to compare the results of the series 
expansion with those of the ICM a similar sampling point distri- 
bution of the projection function (data) was used for both cases 
even though the series expansion method requires neither uniform 
spacing nor as large a number of projections as the ICM. 

When one uses the program SEEX three parameters in addition 
to those used in the ICM program must be specified. They are 
EPSL, N' and M'. Each is defined in Sec. 8.1. In all cases 
reported here, EPSL ■ 10”^. The parameters N* and M' # define the 
number of Legendre polynomials and harmonics in the series repre- 
sentation of the projection function. One series of computations 
was made with N'» M' - 6 and the other with N* * M* ■ 12. Table 
4 presents the maximum and average errors of the SEEXM (series 
expansion method) for all computed reconstructions. The investi- 
gated test function (TOF), opaque obstacle (OB) and the number of 
sampling point M along each diameter over the projection function 
(data) are indicated in the first three columns of Table 4. The 
fourth and fifth columns present the maximum and average recon- 
struction errors for N ' ■ M ' * 6 and N' ■ M' « 12, respectively. 
The last column is identical with the fourth column of the Table 
2 and is repeated only for more convenient comparison of the 
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results of the SEEXM and those of the ICM. 

For N* ■ M' * 12 the SEEXM provides us with nearly zero 
average as well as maximum error in some cases# however# in other 
cases this method completey fails# yielding a 100% maximum error 
and an average error of 15%. For N' * M' * 6 the results of the 
SEEXM were worse or roughly the same as those of the ICM. From 
this very limited investigation we can only conclude that* 

1) For a small number of approximating Fourier harmonics 
and Legendre polynomials the error of the reconstruction is not 
as strongly dependent on the structure of the object field as in 
the case of higher values of N* and M'. 

2) For lower N* and M' the SEEXM provides reconstruction 
which at best are as good as those using the ICM. 

3) The chance that the SEEXM with a higher N' and M‘ will 
give better results than the ICM is good if the investigated 
object function is a simple function only with one maximum or 
minimum and if the area occupied by an opaque obstacle is smaller 
than 50% of the object region. 

Finally# we note that it would be interesting to investigate 
how the reconstruction error is affected by an unequal 
distribution of sampling points of the projection function (data). 

10. The Reconstruction of Data from Aerodynamic Experiments by SEEX 

The program SEEX was used to reconstruct aerodynamic density 
fields from NASA-supplied inter ferograms in the same manner as 
described in Sec. 7 for the ICOM. Reconstructions were computed 
in the three cross-sections z - 0, 1.27 and 2.54cm for flows at 
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TABLE 4 

Comparison of the Maximum and Average Percent 
Errors of the SEEXM and the ICM 


TOF 

OB 

M 

N ' «M ' -6 

N'«M'«12 

ICM 

NOlb 

C3 

21 

8.7 

0.05 

1.5 




1.8 

0.01 

0.2 


C6 

21 

5.0 

100.0 

5.7 




1.1 

19.5 

0.9 


TS 

20 

11.5 

0.08 

14.4 




1.8 

0.01 

1.6 


TB 

20 

13. 2 

0.06 

22.2 




1.5 

0.01 

2.4 

NOld 

C3 

21 

7.8 

0.05 

6.1 




1.8 

0.01 

0.8 


C6 

21 

8.5 

100.0 

11.0 




1.9 

15.9 

1.9 


TS 

20 

8.8 

0.07 

14.9 




1.8 

0.01 

1.8 


TB 

20 

9. 1 

0.06 

19.2 




1.9 

0.01 

2.4 

NO 2 

C3 

27 

47. 5 

25. 1 

19. 5 




10.7 

6.1 

1.8 


C6 

27 

60.0 

100.0 

27. 5 




11.0 

16.9 

3.7 


TS 

20 

43.2 

52.4 

16. 1 




10.3 

8.8 

2.9 


TB 

20 

50.1 

54.9 

22. 6 




in 

« 

c. 

H 

8.4 

3.8 


167 

Bfatfti.1ti.aiii.riif. t i irV-iiT— 


iaeaagfcSssaKasgss^^Mi 


both M ■ 0.6 and M * 0.8. The aerodynamic test object was set at 
a 5.5° angle of attack. {See Sec. 7.) 

Input data for SEEX have been defined in exactly the same 
way as for the reconstructions of the ICM. This means that the 
same interpolations and extrapolations of the measured data and 
the same sampling points were used for input data for the program 
SEEX as for the programs INTEP or ICOM. The parameter EPSL was 
given the value 10“ 4 and the number N' of Legendre polynomials 
and number M' of Fourier harmonics was 12 in all the computed 
cases. The output specifications were the same as for the 
reconstructions by the ICM. Each reconstruction was computed 
with the optional output given by MPLOT > 0 as well as for MPLOT 
0. All the inputs and outputs of the program SEEX are presented 
in the Appendix pp. A107 to A154. 

Density distributions reconstructed using the series 
expansion method are presented graphically in Figs. 13 and 14 for 
Mach numbers M = 0.6 and M * 0.8, respectively. The dashed lines 
are used to display the reconstructed density change Ap(R, <j>). No 
results reconstructed by the SEEXM are presented for the cross- 
sections z = 0 and 1.27cm at Mach number M * 0.8. For these two 
cases the SEEXM completely failed. The reconstructed density 
difference was nearly zero throughout the region. 

The following summarizes a comparison of the reconstructions 
by the series expansion method (SEEXM) and the iterative 
convolution method (ICM): 

1) The reconstructions from complete projection data at the 
tangential plane z * 0 are essentially identical by the two 
methods . 
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2) The discrepancy between results of the two methods at 
other cross-sections is greater than one would expect in view of 
the accuracy found with analytical test fields. 

3) The fields reconstructed by the ICM fluctuate more with 
the angular variable than those reconstructed by the SEEXM. 

4) For fluctuations along the radial coordinate R the 
converse is true. 

The difference in results obtained by the two methods was 
substantially greater than would be expected on the basis of 
reconstruction of analytical test fields. This is probably 
because the real data were noisier and less consistent than the 
simulated data. The interferometric data from the experiment 
were smoothed in the radial direction by visual interpolation, 
but were not smoothed in the azimuthal direction. This led to a 
reconstructed density distribution which was too corrugated 
azimuthal ly when the iterative convolution method was used. The 
series expansion method resulted in more smoothing in this direc- 
tion. Although the SEEXM results in smoothing of the projection 
functions through solving an overdetermined system of equations 
by a least squares method, false fluctuations of the field occur 
in the nodial direction due to non-optimal selections of N' and 
M* . 

Finally, and perhaps most importantly, we note that the 
inter feroqrams used in this study were not of sufficient quality 
to permit accurate reconstruction of the density distribution. 
Because the fringe shifts were not large, the substantial amount 
of background i ;ise forced us to use far too much visual 
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extrapolation as well as ad hoc procedures to estimate fringe 
order* Improvements in experimental technique and fringe 
analysis will be required to obtain accurate reconstructions. 

11 . Conclusions 

We have developed a new type of tomographic code which is 
intended to produce good reconstructions when a test model 
obscures part of the field of view when an inter ferogram of an 
aerodynamic flow field is made. The code combines the features 
of convolution-type algorithms with those of interative algebraic 
techniques. In particular, reconstruction is performed during a 
sequence of iterations between the object field space and the 
Radon transform space (projection space). This iterative scheme 
produces an artificial density field within a convex region which 
bounds the obscuring object. Everywhere outside this obscured 
region the projections, or Radon transform, are represented by 

t 

the measured data from the experiment. The algorithm is 
relatively fast and storage efficient, as is desirable for 
handling large amounts of experimental data. Computer codes 
based on this algorithm have been run on a large microcomputer 
(LSI 11/23) as well as on a mainframe computer. 

This algorithm was studied empirically by using it to 
reconstruct "data" which were line integrals of analytical 
"object fields." Reconstruction errors were determined by 
comparing the object fields reconstructed by the ICM with exact 
values. The effects of the following parameters and factors were 
studied in this manner: 
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i. M, the number of sampled projection values in each 
viewing direction. 

ii. N, the number of views ( inter ferograms ) . 

iii. Technique for initializing data in the region of 
blocked projections. 

iv. Weighting function to be used in the convolution 
algorithm. 

v. Numerical integration method to be used for 
computing projections. 

vi. Criterion for terminating the iterative recon- 
struction procedure. 

The remaining parameters and factors were studied by a 
series of numerical experiments in which two "object fields" were 
reconstructed. One is a field having four local maxima, and the 
other is a field having one maximum located close to its boundary. 
Both fields were zero outside a circle of radius r Q . Four opaque 
objects (i.e. cross-section of aerodynamic models) were used in 
the study: circles of radii r/r Q ■ 0.3 and r/r Q * 0.6, and two 

triangles. Three possible techniques for initializing data (iii) 
were considered. Some analytical studies were used to provide 
guidance for selecting the ratio N/M. Criteria for termination 
(vi) were developed using two figures of merit: the RMS error in 

the reconstruction on the boundary r ■ r Q , where the object field 
must vanish and the RMS error in calculated projections relative 
to their "measured" values over the accessible domain. 

The results of these studies are displayed in a sequence of 
65 figures which are included in this report. Maximum 
reconstruction errors ranged from 1.5% to 47% in these tests and 


171 


average reconstruction errors ranged from 0.2% to 5.3%, depending 
on the complexity of the field, size of the opaque object, etc. 

From this study we developed good guidelines for the 
selection of M and N, the best technique for initializing daca in 
the region of blocked projections, and a fairly reliable 
criterion for terminating the iterative process. 

For comparison, and to develop an altnerate reconstruction 
code applicable to analysis of inter ferograms with blocked 
regions, we wrote a program based on earlier work by Cha in this 
laboratory. This code implicitly fills in the missing 
projections by using a series expansion that meets 
all mathematical criteria for being a Radon transform, in 
particular the k th moment of the projections with respect to the 
radial coordinate with an azimuthal harmonic weighting function 
vanishes if k is less than m, in order of the harmonic weighting 
function. 

When the object field had a simple structure and the opaque 
object occupied less than 50% of the object field area, the 
series expansion gave results comparable with, and in some cases 
more accurate than, the iterative convolution method. When the 
object field was more complicated and/or the opaque object was 
large, the iterative convolution method gave substantially more 
accurate reconstructions. 

Both the iterative convolution and series expansion methods 
were used to reconstruct aerodynamic density fields from 
sequences of holographic inter ferograms of the flow about a 
cylindrical test object with a hemispherical nose. The test 
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object had an angle of attack of 5.5° in an airflow at Mach 
numbers M - 0.6 and M ■ 0.8. The inter ferograms were recorded 
with viewing directions every 10° over a 180° range. It appeared 
that the iterative convolution method gave better 
reconstructions. The series-expansion reconstructions tended to 
be overly corrugated in the radial directions in cross-sectional 
planes that included the test object. The iterative convolution 
reconstructions tended to be corrugated in the nodial direction. 

It should be noted that the inter ferograms contained so much 
extraneous fringe structure (experimental measurement error) that 
rather ad hoc procedures had to be used to assign fringe orders, 
and visual extrapolation had to be used to provide much of the 
data for the computer codes. Improvement in experimental 
techniques and/or fringe processing will be required if accurate 
tomographic reconstructions are to be made. 

Finally, some simple experiments were conducted to show that 
data can be obtained by a holographic technique in which 
interferograms are formed using laser light diffusely reflected 
by the opaque test object. This work was reported separately 
[17,18]. This technique potentially could be used to provide the 
equivalent of a full set of data for tomographic reconstruction 
even though an opaque test object is present. Development of this 
method would require further study of the behavior of tomography 
algorithms when the field includes a region of zero density 
change bounded by a discontinuity. 
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